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Abstract 

Improvements  in  GPS  motivate  attention  to  some  small  relativistic  effects 
which  have  long  been  known,  but  have  previously  been  too  small  to  be  explicitly 
considered.  For  SV  clocks,  these  include  frequency  changes  due  to  orbit  adjust¬ 
ments,  and  effects  due  to  the  earth’s  oblateness.  For  example  between  25  July 
and  10  October  2000,  SV43  occupied  a  transfer  orbit  while  it  was  moved  from  slot 
F5  to  slot  F3.  The  fractional  frequency  change  associated  with  a  change  da  in 
the  semi-major  axis  a  (in  meters)  can  be  estimated  as  9.429  x  10-18da.  This  yields 
a  prediction  of  -1.77  x  10-13  for  the  fractional  frequency  change  of  the  SVfS  clock 
which  occurred  25  July  2000.  This  relativistic  effect  has  been  pointed  out  and 
measured  by  Epstein,  Fine,  and  Stoll  [4].  On  October  10,  2000  the  fractional  fre¬ 
quency  change  should  have  been  +1.75  x  10-13.  Also,  the  earth’s  oblateness  causes 
a  periodic  fractional  frequency  shift  with  period  of  almost  6  hours  and  amplitude 
0.695  x  10-u.  These  effects  will  be  discussed  with  the  help  of  Lagrange’s  planetary 
perturbation-  equations. 


INTRODUCTION 

The  importance  of  relativistic  contributions  to  atomic  clock  frequency  shifts  in  the 
Global  Positioning  System  (GPS)  has  been  recognized  from  the  early  design  stages  of 
the  GPS.  Five  distinct  relativistic  effects  are  incorporated  into  the  System  Specification 
Document,  ICD-GPS-200  [1].  These  are:  the  effect  of  the  earth’s  mass  on  gravitational 
frequency  shifts  of  atomic  reference  clocks  fixed  on  the  earth’s  surface  relative  to  clocks 
at  infinity;  the  effect  of  earth’s  oblate  mass  distribution  on  gravitational  frequency 
shifts  of  atomic  clocks  fixed  on  the  earth’s  surface;  second-order  Doppler  shifts  of 
clocks  fixed  on  the  earth’s  surface  due  to  the  earth  rotation;  gravitational  frequency 
shifts  of  clocks  in  GPS  satellites  due  to  the  earth’s  mass;  and  second-order  Doppler 
shifts  of  clocks  in  GPS  satellites  due  to  their  motion  through  an  Earth-Centered  Inertial 
(ECI)  Frame.  The  combination  of  second-order  Doppler  and  gravitational  frequency 
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shifts  for  a  clock  in  a  GPS  satellite  leads  to  the  following  expression  for  the  fractional 
frequency  shift  of  a  satellite  clock  relative  to  a  reference  clock  fixed  on  the  earth’s 
geoid  [2]: 


Sf  =  _lt£  _  GM  _  $o  m 

/  2  c2  rc 2  c2 

where  v  is  the  satellite  speed  in  a  local  ECI  reference  frame,  GM  is  the  product  of  the 
Newtonian  gravitational  constant  G  and  the  earth’s  mass  M,  c  is  the  defined  speed 
of  light,  and  4>0  is  the  effective  gravitational  potential  on  the  earth’s  rotating  geoid. 
The  term  $o  includes  contributions  from  both  monopole  and  quadrupole  moments  of 
the  earth’s  mass  distribution,  and  the  effective  centripetal  potential  in  an  earth-fixed 
reference  frame  such  as  WGS-84,  due  to  the  earth’s  rotation.  For  reference  we  quote 
here  the  value  for  4>0[2]: 


^  [1  +  J2/2]  -  -Ln \a\  =  -6.9692842  x  10"10  ,  (2) 

c2  aic2  2c2 

where  oj  is  the  earth’s  equatorial  radius,  J2  is  the  quadrupole  moment  coefficient  of 
the  earth,  and  Qe  is  the  angular  rotational  speed  of  the  earth. 

If  the  GPS  satellite  orbit  can  be  approximated  by  a  Keplerian  orbit  of  semi-major  axis 
a,  then  at  an  instant  when  the  distance  of  the  clock  from  the  earth’s  center  of  mass  is 
r,  this  leads  to  the  following  expression  for  the  fraction  frequency  shift  of  Eq.  (1): 

A/  _  3 GM  _  £0  2GM 

f  2ac2  c2  c2 

Eq.  (3)  is  derived  by  making  use  of  the  conservation  of  total  energy  (per  unit  mass) 
of  the  satellite,  which  leads  to  an  expression  for  v2  in  terms  of  GM/r  and  GM/a  which 
can  be  substituted  into  Eq.  (1): 

1  2  GM  GM 

2  r  2a  w 

The  first  two  terms  in  Eq.  (3)  give  rise  to  the  “factory  frequency  offset,”  which  is 
applied  to  GPS  clocks  before  launch  in  order  to  make  them  beat  at  a  rate  equal  to 
that  of  reference  clocks  on  the  earth’s  surface.  The  last  term  in  Eq.  (3)  is  very  small 
when  the  orbit  eccentricity  e  is  small;  when  integrated  over  time,  these  terms  give 
rise  to  the  so-called  ‘e  sin  E”  effect  or  “eccentricity  effect.”  In  most  of  the  following 
discussion,  we  shall  assume  that  eccentricity  is  very  small. 

Clearly,  from  Eq.  (3),  if  the  semi-major  axis  should  change  by  an  amount  Sa  due  to 
an  orbit  adjustment,  the  satellite  clock  will  experience  a  fractional  frequency  change 

6f  ,  3 GM6a 

f  ~  +  2 c2a2  '  (  } 

The  factor  3/2  in  this  expression  arises  from  the  combined  effect  of  second-order 
Doppler  and  gravitational  frequency  shifts.  If  the  semi-major  axis  increases,  the  satel¬ 
lite  will  be  higher  in  the  earth’s  gravitational  potential  and  will  be  gravitationally 
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blue-shifted  more,  while  at  the  same  time  the  satellite  velocity  will  be  reduced,  reduc¬ 
ing  the  size  of  the  second-order  Doppler  shift  (which  is  generally  a  red  shift).  The  net 
effect  would  make  a  positive  contribution  to  the  fractional  frequency  shift. 

It  has  long  been  known  that  orbit  adjustments  are  associated  with  satellite  clock 
frequency  shifts  [3],  but  nothing  has  been  documented  and  no  reliable  measurements 
of  such  shifts  have  previously  been  made.  Recently  Marvin  Epstein,  Joseph  Fine, 
and  Eric  Stoll  [4]  of  ITT  have  evaluated  the  frequency  shift  of  SV43  arising  from  a 
trajectory  change  applied  to  this  satellite  on  25  July  2000.  The  purpose  of  this  orbit 
adjustment  was  to  move  the  satellite  from  slot  F5  to  slot  F3.  A  drift  orbit  extending 
from  25  July  2000  to  10  October  2000  was  used  to  accomplish  this  move  [6].  Epstein, 
Fine,  and  Stoll  have  reported  that  associated  with  the  “delta-v  burn”  on  25  July  2000 
there  was  a  frequency  shift  of  the  rubidium  clock  on  board  SV43  of  amount 

=  — 1.85  x  10“ 13  (measured) .  (6) 

Epstein  et  al. [4]  suggested  that  the  above  frequency  shift  is  relativistic  in  origin,  and 
used  NIMA  precise  ephemerides  to  estimate  the  frequency  shift  arising  from  second- 
order  Doppler  and  gravitational  potential  differences.  They  calculated  separately  the 
second-order  Doppler  and  gravitational  frequency  shifts  due  to  the  orbit  change.  The 
NIMA  precise  ephemerides  are  expressed  in  the  WGS-84  coordinate  frame,  which  is 
earth-fixed.  If  used  without  removing  the  underlying  earth  rotation,  the  velocity  would 
be  erroneous.  They,  therefore,  transformed  the  NIMA  precise  ephemerides  to  an  earth- 
centered  inertial  frame  by  accounting  for  a  (uniform)  earth  rotation  rate. 

The  semi-major  axes  before  and  after  the  orbit  change  were  calculated  by  taking  the 
average  of  the  maximum  and  minimum  radial  distances.  Speeds  were  calculated  using 
a  Keplerian  orbit  model.  They  [4]  arrived  at  the  following  numerical  values  of  position 
and  velocity: 


07/22/00  :  a  =  2.656139556  x  107  m.;  v  =  3.873947951  x  103  m/s.  (7) 

07/30/00  :  a  =  2.654267359  x  107  m.;  v  =  3.875239113  x  103  m/s.  (8) 

Since  the  semi-major  axis  decreased,  the  frequency  shift  should  be  negative.  The 
prediction  [4]  made  for  the  frequency  shift,  which  is  based  on  Eq.  (1),  was  then: 

^  =  -1.734  x  10"13  (9) 

which  is  to  be  compared  with  the  measured  value,  Eq.  (6).  There  is,  thus,  fairly 
compelling  evidence  that  the  observed  frequency  shift  is  indeed  a  relativistic  effect. 


LAGRANGE  PERTURBATION  THEORY 


Perturbations  of  GPS  orbits  due  to  the  earth’s  quadrupole  mass  distribution  are  a 
significant  fraction  of  the  change  in  semi-major  axis  associated  with  the  orbit  change 
discussed  above.  This  raises  a  question  whether  it  is  sufficiently  accurate  to  use  a 
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Keplerian  orbit  to  describe  GPS  satellite  orbits,  and  estimate  the  semi-major  axis 
change  as  though  the  orbit  were  Keplerian.  The  purpose  of  the  present  calculation  is 
to  investigate  this  question.  Also,  it  is  of  interest  to  investigate  directly  the  effect  of 
the  earth’s  quadrupole  moment  on  the  frequency  of  GPS  satellite  clocks.  Previously, 
such  effects  have  been  neglected.  However,  the  effect  may  be  worth  considering  as 
GPS  clock  performance  continues  to  improve. 

To  see  how  large  such  quadrupole  effects  may  be,  the  Appendix  [5]  quotes  the  perturba¬ 
tions  on  the  semi-major  axis  a  of  a  Keplerian  satellite  arising  from  earth’s  quadrupole 
moment.  For  the  semi-major  axis,  if  the  eccentricity  is  very  small  the  dominant  contri¬ 
bution  has  a  period  twice  the  orbital  period  and  has  amplitude  3J2a2  sin2*o/(2ao)  «  1658  m. 
The  following  values  for  the  constants  are  used  here: 

J2  =  1.0826267  x  1(T3;  GM  =  3.986004415  x  1014m3/  sec2;  m  =  6.3781363  x  106  m;  QE  =  7.291151467  x 
KT5  sec-1;  a0  =  2.65641046  x  107  m, 

where  a0  and  ai  are  the  orbit  semi-major  axis  and  the  earth’s  equatorial  radius,  re¬ 
spectively,  and  He  is  the  earth’s  rotational  angular  velocity. 

The  oscillation  in  the  semi-major  axis  would  significantly  affect  calculations  of  the 
semi-major  axis  at  any  particular  time.  This  suggests  that  Eq.  (4)  needs  to  be  exam¬ 
ined  carefully  in  light  of  the  periodic  perturbations  on  the  semi-major  axis.  Therefore 
in  this  paper  we  develop  an  approximate  description  of  a  satellite  orbit,  of  small  ec¬ 
centricity,  taking  into  account  the  earth’s  quadrupole  moment  to  first  order.  Terms  of 
order  J2xe  will  be  neglected.  This  problem  is  non-trivial  because  the  perturbations 
themselves  (see  for  example,  the  equations  for  mean  anomaly  and  altitude  of  perigee) 
have  factors  1/e  which  blow  up  as  the  eccentricity  approaches  zero.  This  problem  is  a 
mathematical  one,  not  a  physical  one.  It  simply  means  that  the  observable  quantities- 
such  as  coordinates  and  velocities-need  to  be  calculated  in  such  a  way  that  finite  values 
are  obtained.  Orbital  elements  which  blow  up  are  unobservable. 


CONSERVATION  OF  ENERGY 


The  gravitational  potential  of  a  satellite  at  position  (x,y,z)  in  equatorial  ECI  coordi¬ 
nates  in  the  model  under  consideration  here  is 


V(x,y,z)  =  - 


GM 

r 


J2a2  '3z2 
r2  2  r2 


(10) 


Since  the  force  is  conservative  in  this  model  (solar  radiation  pressure,  thrust,  etc.  are 
not  considered),  the  kinetic  plus  potential  energy  is  conserved.  Let  e  be  the  energy  per 
unit  mass  of  an  orbiting  mass  point.  Then 


v2  x  v2  GM  . 

e  =  constant  =  —  +  V(x,  y,  z)  =  — - f-  V  ( x,y,z ) 


(11) 


where  V'(x,y,z)  is  the  perturbing  potential  due  to  the  earth’s  quadrupole  potential. 

It  is  shown  in  textbooks  [6]  that,  with  the  help  of  Lagrange’s  planetary  perturbation 
theory,  the  conservation  of  energy  condition  can  be  put  in  the  form 
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(12) 


Kjrivi  .  x 

e=-^-+V(x,y,z) 

where  a  is  the  perturbed  (osculating)  semi-major  axis,  given  by  Eq.  (A. 6)  in  the 
Appendix.  In  other  words,  for  the  perturbed  orbit, 

v2  GM  _  GM 

2  r  ~  2a  1  ' 

On  the  other  hand,  the  net  fractional  frequency  shift  relative  to  a  clock  at  rest  at 
infinity  is  determined  by  the  second-order  Doppler  shift  (a  redshift)  and  a  gravitational 
redshift.  The  total  relativistic  fractional  frequency  shift  is 


A / 
/ 


GM 

r 


+  V'(x,y,z) 


(14) 


The  conservation  of  energy  condition  can  be  used  to  express  the  second-order  Doppler 
shift  in  terms  of  the  potential.  Since  in  this  paper  we  are  interested  in  fractional 
frequency  changes  caused  by  changing  the  orbit,  it  will  make  no  difference  if  the  cal¬ 
culations  use  a  clock  at  rest  at  infinity  as  a  reference  rather  than  a  clock  at  rest  on  the 
earth’s  surface.  The  reference  potential  cancels  out  to  the  required  order  of  accuracy. 

Therefore,  from  perturbation  theory  we  need  expressions  for  the  square  of  the  velocity, 
for  the  radius  r,  and  for  the  perturbing  potential.  We  now  proceed  to  derive  these 
expressions. 


PERTURBATION  EQUATIONS 

First  we  list  some  facts  about  an  unperturbed  Keplerian  orbit.  The  eccentric  anomaly 
E  is  to  be  calculated  by  solving  the  equation 


E  —  e  sin  E  =  M  =  n0(t  —  to)  (15) 

where  M  is  the  “mean  anomaly”  and  t0  is  the  time  of  passage  past  perigee,  and 


n0  ~  \J GAf/a 3  .  (16) 

Then  the  perturbed  radial  distance  r  and  true  anomaly  f  of  the  satellite  are  obtained 
from 


r  =  a(l  —  ecosE) 


(17) 


cos  /  = 


cos  E  —  e 
1  —  e  cos  E  ’ 


sin  /  =  x/T 


sin  E 


1  —  e  cos  E 


(18) 


The  observable  £,y,z-coordinates  of  the  satellite  are  then  calculated  from  the  following 
equations: 
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x  =  r(cosfJcos(/  +  a;)  —  cosisinfisin(/  4-  a;)) 


(19) 


y  =  r(sinficos(/  +  w)  +  cosicosfisin(/  +  a>)). 


(20) 


z  =  r(sinisin(/  +  w))  (21) 

where  a;  is  the  angle  of  the  ascending  line  of  nodes,  i  is  the  inclination,  and  w  is  the 
altitude  of  perigee. 

By  differentiation  with  respect  to  time,  or  by  using  the  conservation  of  energy  equation, 
one  obtains  the  following  expression  for  the  square  of  the  velocity: 

v 2  _  GM  1  +  ecos  E  .  . 

2  2a  1  —  ecos  E 

In  these  expressions  v2  and  r-1  are  observable  quantities.  The  combination  ecos  E,  where 
E  is  the  eccentric  anomaly,  occurs  in  both  of  these  expressions.  To  derive  expressions 
for  v2  and  r-1  in  the  perturbed  orbits,  expressions  for  the  perturbed  elements  a,  e,  E 
are  to  be  substituted  into  the  righthand  sides.  Therefore,  we  need  the  combination 
ecos E  in  the  limit  of  small  eccentricity. 


Calculation  of  Perturbed  Eccentricity 

To  leading  order,  from  the  Appendix  we  have  for  the  perturbed  eccentricity  the  fol¬ 
lowing  expression: 


e=Ke  + 


3J2aj 
2  Oq 


2  r 


3  \  i  7 

1  -  -  sin2  i0  )  cos  /  +  -  sin2  i0  cos(2w0  +  /)  +  —  sin2  i0  cos(2w0  +  3 /) 


(23) 


Calculation  of  Perturbed  Eccentric  Anomaly 

The  eccentric  anomaly  is  calculated  from  the  equation 


E  =  M  +  e  sin  E  (24) 

with  perturbed  values  for  M  and  e.  Expanding  to  first  order  in  e  gives  the  following 
expression  for  cos  E: 


cos  E  -  cos  M  —  e  sin  M  sin  E 


(25) 


and  multiplying  by  e, 


e  cos  E  =  e  cos  M  —  e2  sin  M  sin  E  «  e  cos  M 


(26) 
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We  shall  neglect  higher  order  terms  in  e.  Referring  now  to  the  perturbed  expression 
for  mean  anomaly  M  in  the  Appendix,  we  write 


M  =  Mo  +  A  M/e0.  (27) 

where  we  indicate  explicitly  the  terms  in  e^1;  that  is,  the  quantity  Mq  contains  all 
terms  which  do  not  blow  up  as  e  — >  0,  and  AM/ e0  contains  all  the  other  terms.  From 
the  Appendix,  we  have  to  leading  order 


AM/e0  = 


3  J2a? 
2  e0al 


1  7 

sin  /  —  -  sin2  io  sm(2w0  +  /)  +  —  sin2  io  sin(2o>o  +  3/) 


and  so  for  very  small  eccentricity, 


(28) 


e  cos  E  —  e  cos  Mo  —  AM  sin  Mo  ■  (29) 

Then  after  accounting  for  contributions  from  the  perturbed  eccentricity  and  the  per¬ 
turbed  mean  anomaly,  after  a  few  lines  of  algebra  we  obtain  for  ecos  E 

e  cos  E  —  e0  cos  Eq  +  ~7£~T~  ^  sin2  sin2  io  cos  2(o>o  +  /).  (30) 

where  the  first  term  is  the  unperturbed  part. 

The  perturbation  is  a  constant,  plus  a  term  with  twice  the  orbital  period. 

Calculation  of  the  Perturbation  in  Semi-major  Axis 

From  the  Appendix,  the  leading  terms  in  the  perturbation  of  the  semi-major  axis  are 


a  =  Ka  + 


3  J2af 
2<xq 


sin2  io  cos  2  ( wq  +  /)• 


(31) 


Calculation  of  the  Perturbation  in  Radius 

We  are  now  in  position  to  compute  the  perturbation  in  the  radius.  From  the  expression 
for  r,  we  have  after  combining  terms 

r  =  a0(l  —  eocosEo)  +  Aa  —  A(ecos  E) 

=  a0(l  -  e0  cos  Eq)  -  -^Q-  ^1  -  |sin2i0^  +  :^-sin2*ocos2(wo  +  /)  •  (32) 
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Calculation  of  the  Perturbation  in  the  Velocity  Squared 

The  above  results,  after  substituting  into  Eq.  (18),  yield  the  expression 


v 

~2 


GM_ 

2ao 


(1  +  2eo  cos  Eq)  + 


3  GMJ2a\ 


1  -  -  sin2  *0  + 


GMJ^aj 
2oq 


sin2  *o  cos  2(<jJq  +  /) .  (33) 


Calculation  of  the  Perturbation  in  GM/r 

The  above  expression  for  the  perturbed  r  yields  the  following  for  the  monopole  contri¬ 
bution  to  the  gravitational  potential: 


GM_ 

r 


GM_ 

ao 


(1  +  e0  cos  Eq)  - 


3  GMJ2af 
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Calculation  of  the  Approximate  Value  of  the  Perturbing  Potential 

Since  the  perturbing  potential  contains  the  small  factor  J2,  to  leading  order  we  may 
substitute  unperturbed  values  for  r  and  2  into  V'(x,y,z)  which  yields  the  expression 
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(35) 


Conservation  of  Energy 

It  is  now  very  easy  to  check  conservation  of  energy.  Adding  kinetic  energy  per  unit 
mass,  to  two  contributions  to  the  potential  energy,  gives 
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v2  GM  GM 
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2  r  2a0 


GAGhaj 
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(36) 


This  verifies  that  the  perturbation  theory  gives  a  constant  energy.  The  extra  term 
in  the  above  equation,  with  J2  in  it,  can  be  neglected.  This  is  because  the  nominal 
inclination  of  GPS  orbits  is  such  that  the  factor  (1  —  3sin2i0/2)  is  essentially  zero. 

Thus  numerical  calculations  of  the  total  energy  per  unit  mass  should  give  us  the  value 
of  ao. 


Calculation  of  Fractional  Frequency  Shift 

The  fractional  frequency  shift  calculation  is  very  similar  to  the  calculation  of  the 
energy,  except  the  the  second-order  Doppler  term  contributes  with  a  negative  sign. 
The  result  is 
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3  GM_  2  GM 

2  a0c2 


O - e0  COS  Eq  — 

c2a0 


TGMJzaj 

2oqC2 


(l-^sin2i0) 


GMJ^a2  sin2  *o 


«oc 2 
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The  first  term,  when  combined  with  the  reference  potential  at  the  earth’s  geoid  gives 
rise  to  the  “factory  frequency  offset.”  The  third  term  can  be  neglected,  as  pointed  out 
above.  The  last  term  has  an  amplitude 


ale 2 


(38) 


which  may  be  large  enough  to  consider  when  calculating  frequency  shifts  produced  by 
orbit  changes.  Therefore,  this  contribution  may  have  to  be  considered  in  the  determi¬ 
nation  of  the  semi-major  axis. 

The  result  suggests  the  following  method  of  computing  the  fractional  frequency  shift: 
Averaging  the  shift  over  one  orbit,  the  periodic  term  will  average  down  to  a  negligible 
value.  The  second  term  is  negligible.  So  if  one  has  a  good  estimate  for  the  nominal 
semi-major  axis  parameter,  the  term  -3GM/2aoc2  gives  the  average  fractional  frequency 
shift.  On  the  other  hand,  the  average  energy  per  unit  mass  is  given  by  e  =  -GM/2a0. 
Therefore,  the  precise  ephemerides,  specified  in  an  ECI  frame,  can  be  used  to  compute 
the  average  value  for  e,  then  the  average  fractional  frequency  shift  will  be 


t=3£/c2' 


(39) 


The  periodic  term  in  Eq.  (33)  is  of  a  form  similar  to  that  which  gives  rise  to  the  ec¬ 
centricity  correction,  which  is  applied  by  GPS  receivers.  Considering  only  the  periodic 
term,  the  additional  time  elapsed  on  the  orbiting  clock  will  be  given  by 


Stj2  =  j  dt 

J  path 


QV/Wsin  >cos(2ut,  +  2„() 
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(40) 


where  to  a  sufficient  approximation  we  have  replaced  the  quantity  2/  in  the  integrand 
by  2 n  =  2y/GM/a$;  n  is  the  approximate  mean  motion  of  GPS  satellites.  Upon  integrat¬ 
ing  and  dropping  the  constant  of  integration  (assuming  as  usual  that  such  constant 
time  offsets  are  lumped  with  other  contributions)  gives  the  periodic  relativistic  effect 
on  the  elapsed  time  of  the  SV  clock  due  to  earth’s  quadrupole  moment: 


GM  J2a\  sin2  i0  , 

StJ2  =  - ^2 - sm(2w0  +  2 nt) . 


(37) 


The  correction  which  should  be  applied  by  the  receiver  is  the  negative  of  this  expression 


c.  ,  .  / GM  J2a\  sin2  i0  .  „  .. 

otj^ (correction)  =  \  j— 3 - ^2 - sin(2a>o  +  2nt) . 


(38) 
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The  phase  of  this  correction  is  zero  when  the  satellite  passes  through  the  earth’s 
equatorial  plane  going  northwards. 

If  not  accounted  for,  this  effect  on  the  SV  clock  time  would  give  rise  to  a  peak-to-peak 
periodic  navigational  error  in  position  of  approximately  2 cxStj2  =  1.07  cm. 

These  effects  were  considered  by  Ashby  and  Spilker  ([2],  pp.  685-686]),  but  in  that  work 
the  effect  of  the  earth’s  quadrupole  moment  on  the  term  GM/r  was  not  considered; 
the  present  paper  supersedes  that  work. 


NUMERICAL  CALCULATIONS 

Precise  ephemerides  were  obtained  for  SV43  from  the  Jet  Propulsion  Laboratories  Web 
site  ftp://sideshow.jplnasa.gov/pub/2000/orbits .  These  are  expressed  in  the  J2000  ECI 
frame.  Computer  code  was  written  to  compute  the  average  value  of  e  for  one  day  and 
thence  the  fractional  frequency  shift  relative  to  infinity  before  and  after  each  orbit 
change.  The  following  results  were  obtained: 


07/22/00  :  a  -  2.65611575  x  107  ±  69  m. 

07/30/00  :  a  =  2.65423597  x  103  ±  188  m. 

10/07/00  :  a  =  2.65418742  x  107  dh  95  m. 

10/12/00  :  a  =  2.65606323  x  107  ±  58  in. 

Therefore,  the  fractional  frequency  change  produced  by  the  orbit  change  of  July  25  is 
calculated  to  be 


-j-  =  -1.77  x  10“13  ,  (36) 

which  agrees  with  the  measured  value  to  within  about  3.3%.  We  predict  that  the 
fractional  frequency  shift  on  October  10,  should  have  been 

=L  = +1.75  x  10-13.  (37) 

This  shift  has  not  yet  been  measured. 

On  9  March  2001,  SV54’s  orbit  was  changed  by  a  delta-u  burn.  Using  the  above 
procedures,  we  can  calculate  the  fractional  frequency  change  produced  in  the  onboard 
clocks.  We  find 


03/07/01  :  a  =  2.65597188  x  107  ±  140  m. 

03/11/01  :  a  =  2.65359261  x  107  ±  108  m. 

Using  Eq.  (5)  yields  the  following  prediction  for  the  fractional  frequency  change  of 
SV54  on  9  March  2001: 
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=  +2.24  x  10-13  ±  0.02  x  10"13 . 

The  quoted  uncertainty  is  due  to  the  combined  uncertainties  from  the  determination 
of  the  energy  per  unit  mass  before  and  after  the  orbit  change. 

CONCLUSIONS 

We  note  that  the  values  of  semi- major  axis  reported  by  Epstein  et  al..  differ  from 
the  values  obtained  by  averaging  as  outlined  above,  by  200-300  m.  This  difference 
arises  because  of  the  different  methods  of  calculation.  In  the  present  calculation,  an 
attempt  was  made  to  account  for  the  effect  of  the  earth’s  quadrupole  moment  on  the 
Kepler ian  orbit.  It  was  not  necessary  to  compute  the  orbit  eccentricity.  Agreement 
with  measurement  of  the  fractional  frequency  shift  was  only  a  few  percent  better  than 
that  obtained  by  differencing  the  maximum  and  minimum  radii. 

This  approximate  treatment  of  the  orbit  makes  no  attempt  to  consider  perturbations 
that  are  non-gravitational  in  nature-e.g.,  solar  radiation  pressure.  The  work  was  an 
investigation  of  the  approximate  effect  of  the  earth’s  quadrupole  moment  on  the  GPS 
satellite  orbits,  for  the  purpose  of  (possibly)  accurate  calculations  of  the  fractional 
frequency  shifts  which  result  from  orbit  changes. 

As  a  general  conclusion,  the  fractional  frequency  shift  can  be  estimated  to  very  good 
accuracy  from  the  expression  for  the  “factory  frequency  offset.” 
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APPENDIX 


This  Appendix  quotes  for  convenience  the  results  of  first-order  Lagrangian  perturba¬ 
tion  theory  for  a  Keplerian  orbit  perturbed  by  a  mass  quadrupole  centered  at  the 
origin  [6],  The  osculating  elements  are:  semi-major  axis  a,  eccentricity  e,  inclination 
i,  longitude  of  the  ascending  node  Q,  altitude  of  perigee  u>,  and  mean  anomaly  M.  La¬ 
grange’s  planetary  perturbation  equations  describe  how  these  six  elements  change  in 
time  under  the  influence  of  a  perturbation.  In  the  present  case,  the  perturbation  is 
due  to  the  earth’s  quadrupole  mass  distribution, 


V'(x,y,z) 


GMjhaj 

1*3 


\3z2 

2r2 


(Al) 


where  GM  is  the  product  of  the  Newtonian  gravitational  constant  and  the  earth’s  mass, 
J2  is  the  earth’s  quadrupole  moment  coefficient,  r  is  the  satellite’s  radial  distance  from 
the  earth’s  center  of  mass,  and  ( x,y,z )  are  the  satellite  coordinates  in  an  earth-fixed 
inertial  reference  frame  with  z-axis  parallel  to  the  earth’s  symmetry  axis. 

When  these  perturbed  orbital  elements  are  given,  the  eccentric  anomaly  E  is  to  be 
calculated  by  solving  the  equation 


E  —  e  sin  E  =  M  ( A2) 

and  then  the  perturbed  radial  distance  r  and  true  anomaly  /  of  the  satellite  are  ob¬ 
tained  from 
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The  observable  x,  y,  z-coordinates  of  the  satellite  are  then  calculated  from  the  following 
equations: 


x  =  r(cos  fi  cos(/  +  u>)  —  cos?,  sin  f2sin(/  +  u>)) 


y  =  r(sin  cos (/  +  w)  +  cos  i  cos  Cl  sin (/  +  a>)). 


(A5) 


z  =  r(sin  i  sin(f  -f-w)) 

The  following  expressions  give  the  perturbed  orbital  elements  correct  to  first  order  in 
the  perturbation. 
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The  perturbed  mean  anomaly  is 


Jl/  —  K M  "t  Tift  — 
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and  e  is  the  conserved  energy  per  unit  mass,  given  by 


GM  v2  T„, 
e= - —  +  ~2+v  (x,y,z) 
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Constants  such  as  Ka,  Ke,  etc.,  are  constants  of  integration  determined  by  the  initial 
conditions. 


Questions  and  Answers 


VICTOR  REINHARDT  (Boeing  Satellite  Systems):  Can  you  put  up  that  next  to  the  last  slide  that  had  the 
number  of  picoseconds  change  due  to  the  varying  term?  I  just  missed  the  slides  of  that.  So  that  was  about 
24  picoseconds. 

NEIL  ASHBY :  Twenty-four  picoseconds. 

REINHARDT:  Okay,  thank  you. 

EDOARDO  DETOMA  (Alenia  Aerospazio):  If  I  remember  well,  at  the  beginning,  the  GPS  inclination  was 
64.3  degrees  or  something  like  that,  not  55.  And  64.3  inclination  is  the  0,1  term  of  the  series  expansion  of 
the  simulation  - 1  believe  it  was  not  incidental;  I  believe  it  was  done  on  purpose. 

ASHBY:  I  think  that  zeroes  out  the  procession  of  the  modal  line. 

DETOMA:  Yes,  but  probably  more  than  that.  There  is  a  paper  which  was  published,  I  believe,  in  ’57  or 
’58  in  which  the  term  for  64.3  is  clearly  pinpointed. 

ASHBY:  Yes.  But  the  banishing  of  this  particular  term  does  not  zero  out  any  secular  terms.  I  don’t  know 
what  the  reason  was.  The  55. 

DETOMA:  At  a  certain  point,  the  constellation  was  reduced  to  18  satellites.  If  you  look  at  the  literature 
around  ’86,  ’87,  you  find  they  change  the  number  of  satellites  from  24  to  18.  And  this  was  probably  due  to 
budgetary  reasons.  And  if  you  reduce  the  number  of  satellites,  you  must  reduce  the  inclination  to  maintain 
the  filling  factor  of  the  constellation. 

DEMETRIOS  MATSAKIS  (U.S.  Naval  Observatory):  I  would  like  to  just  confirm  that  a  little  bit.  It  is  a 
rumor,  but  I  was  told  by  Colonel  Armour  years  ago  that  the  reason  for  the  inclination  was  budget.  It  cost 
more  to  have  a  higher  inclination. 

TOM  McCASKDLL  (U.S.  Naval  Research  Laboratory):  Going  way  back  to  around  1968-1969,  whenever 
Roger  Easton,  Jim  Buisson,  Don  Lynch,  and  myself  did  the  constellation  study,  we  looked  at  about  107 
different  constellation  configurations.  And  we  calculated  a  lot  of  quantities,  and  it  turned  out  our 
calculation  and  that  you  could  minimize  the  dilution  precision  with  an  inclination  of  53  degrees.  And  later 
on,  with  the  Block  I’s,  it  went  up  to  63.  And  there  was  a  study  done  by  Aerospace  and  the  Air  Force,  and 
they  finally  settled  on  55.  And  there  could  have  been  some  launch  considerations;  however,  if  you  look  at 
zero-degree  inclination  for  the  constellation,  there  is  a  mathematical  singularity,  and  if  you  go  up  to  90 
degrees,  like  Transit,  you  have  too  many  satellites  at  the  pole.  So  it  is  really  a  balancing  operation  that  gives 
you  the  best  dilution  of  precision  on  a  worldwide  basis. 
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